1 Reading data

#read data obtained in "01_data_processing.rmd":
load("processed_data.rda")
load("cue_rates_obtained.rda")
sampledata<-sampledata[sampledata$rsecs>=0&!is.na(sampledata$clicks),]
sampledata$time11<-as.numeric(sampledata$rsecs)
originaldata$depth<-originaldata$Depth
originaldata$time11<-as.numeric(originaldata$sst)
onedata$depth<-onedata$Depth
onedata$time11<-as.numeric(onedata$sst)

2 GAM Models

2.1 GAM number of clicks per second as a function of depth by whale and as a funtion of time by whale

gamwhales1<-gam(clicks~s(depth,k=7,by=as.factor(whale))+s(time11,k=1,by=as.factor(whale)),data=sampledata,family=poisson(link="log"))
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible

2.1.1 Summary of the model

summary(gamwhales1)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale)) + s(time11, k = 1, 
##     by = as.factor(whale))
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.39154    0.02194   63.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df Chi.sq  p-value    
## s(depth):as.factor(whale)Freya   4.325  4.940 17.035 0.003860 ** 
## s(depth):as.factor(whale)Thora   5.644  5.940 21.738 0.000932 ***
## s(depth):as.factor(whale)Mára    2.210  2.671  1.631 0.598642    
## s(depth):as.factor(whale)Frida   1.000  1.000 19.735 9.49e-06 ***
## s(depth):as.factor(whale)Eistla  2.349  2.919 16.056 0.000969 ***
## s(depth):as.factor(whale)Balder  4.242  4.999 33.387 3.60e-06 ***
## s(depth):as.factor(whale)Jonas   4.750  5.447 29.933 4.59e-05 ***
## s(depth):as.factor(whale)Mutti   1.000  1.000 21.955 3.52e-06 ***
## s(time11):as.factor(whale)Freya  1.729  1.925 13.614 0.005081 ** 
## s(time11):as.factor(whale)Thora  1.730  1.926  4.325 0.075160 .  
## s(time11):as.factor(whale)Mára   1.000  1.000  0.399 0.527811    
## s(time11):as.factor(whale)Frida  1.396  1.634 24.510 9.76e-05 ***
## s(time11):as.factor(whale)Eistla 1.000  1.000  2.812 0.093580 .  
## s(time11):as.factor(whale)Balder 1.000  1.000 17.759 2.56e-05 ***
## s(time11):as.factor(whale)Jonas  1.001  1.002  0.159 0.690277    
## s(time11):as.factor(whale)Mutti  1.000  1.000  0.369 0.543610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0464   Deviance explained = 6.19%
## UBRE = 1.8611  Scale est. = 1         n = 1203

2.1.2 AIC

AIC(gamwhales1)
## [1] 6732.858

2.1.3 Predicting influence of the different expl. variables

dd1<-sampledata%>%
  group_by(whale)%>%
  summarise(depth=mean(depth,na.rm=T))# creating a data set to predict the influence of time


dd2<-sampledata%>%
  group_by(whale)%>%
  summarise(time11=mean(time11,na.rm=T))# creating a data set to predict the influence of depth

originaldata11<-data.frame(depth=sampledata$depth,whale=sampledata$whale)# creating a data set to predict the influence of depth

originaldata11<-full_join(originaldata11,dd2)# joining both infos to create a data set to predict the influence of depth
## Joining, by = "whale"
originaldata111<-data.frame(time11=sampledata$time11,whale=sampledata$whale)# creating a data set to predict the influence of time

originaldata111<-full_join(originaldata111,dd1)# joining both infos to create a data set to predict the influence of time
## Joining, by = "whale"
datat1<-data.frame(pclick=predict(gamwhales1,newdata = originaldata11,type="response"),depth=originaldata11$depth,whale=originaldata11$whale)# predicting the influence of depth in the number of clicks produced

confintervals<-data.frame(confintervals=predict(gamwhales1,newdata=(originaldata11),type="response",se.fit=T))# get the corresponding CI

datat1$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat1$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit

datat1$variable<-"Depth (m)"# which variable corresponds

datat11<-data.frame(pclick=predict(gamwhales1,newdata = originaldata111,type="response"),time11=originaldata111$time11,whale=originaldata111$whale)# predicting the influence of time in the number of clicks produced

confintervals<-data.frame(confintervals=predict(gamwhales1,newdata=(originaldata111),type="response",se.fit=T))# get the corresponding CI

datat11$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat11$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit

datat11$variable<-"Time"# which variable corresponds


datat<-full_join(datat1,datat11)#Joinning both infos about the variables 
## Joining, by = c("pclick", "whale", "bigfit", "smallfit", "variable")
datat$value<-ifelse(datat$variable=="Depth (m)",datat$depth,datat$time11)# turn the variable depth into the real depth

datat$value<-ifelse(datat$variable=="Time",datat$value/(60*60*24),datat$value)# turn the variable time into the real time in days

datat$variable1<-ifelse(datat$variable=="Time","Time since regular echolocation clicks were resumed (days)",datat$variable)# new labels


datat$whale<-factor(datat$whale,levels=c("Freya","Thora","Mára","Frida","Eistla","Balder","Jonas","Mutti"))# correct order of the whale  id

2.1.4 Plotting model results/ Figure 5

ggplot(data=datat,aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+facet_wrap(~variable1,scales="free_x", labeller = labeller(variable1 = label_wrap_gen(32)))+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

p1<-ggplot(data=subset(datat,variable1=="Time since regular echolocation clicks were resumed (days)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("")+xlab("Time since regular echolocation clicks \nwere resumed (days)")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(axis.text.y=element_blank(),
      #  axis.ticks.y=element_blank(),
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
       # axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))+ylim(2,12)

p2<-ggplot(data=subset(datat,variable1=="Depth (m)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))+ylim(2,12)

ggpubr::ggarrange(p2,p1,ncol=2,common.legend = T,widths = c(0.5,0.4), legend = "right")

2.1.5 Predicting number of clicks for the clicking periods from the original dataset

onedata$prediction1<-predict(gamwhales1,newdata=(onedata),type="response")# get the prediction from the model

2.1.6 Plot with the predicted number of clicks for the clicking periods from the original dataset / Figure 6

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction1),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction1),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

2.2 GAM number of clicks per second as a function of depth by whale and as a funtion of time by whale

gamwhales11<-gam(clicks~s(depth,k=7)+s(time11,k=1),data=sampledata,family=poisson(link="log"))
## Warning in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots): basis dimension, k, increased to minimum possible

2.2.1 Summary of the model

summary(gamwhales1)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale)) + s(time11, k = 1, 
##     by = as.factor(whale))
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.39154    0.02194   63.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                    edf Ref.df Chi.sq  p-value    
## s(depth):as.factor(whale)Freya   4.325  4.940 17.035 0.003860 ** 
## s(depth):as.factor(whale)Thora   5.644  5.940 21.738 0.000932 ***
## s(depth):as.factor(whale)Mára    2.210  2.671  1.631 0.598642    
## s(depth):as.factor(whale)Frida   1.000  1.000 19.735 9.49e-06 ***
## s(depth):as.factor(whale)Eistla  2.349  2.919 16.056 0.000969 ***
## s(depth):as.factor(whale)Balder  4.242  4.999 33.387 3.60e-06 ***
## s(depth):as.factor(whale)Jonas   4.750  5.447 29.933 4.59e-05 ***
## s(depth):as.factor(whale)Mutti   1.000  1.000 21.955 3.52e-06 ***
## s(time11):as.factor(whale)Freya  1.729  1.925 13.614 0.005081 ** 
## s(time11):as.factor(whale)Thora  1.730  1.926  4.325 0.075160 .  
## s(time11):as.factor(whale)Mára   1.000  1.000  0.399 0.527811    
## s(time11):as.factor(whale)Frida  1.396  1.634 24.510 9.76e-05 ***
## s(time11):as.factor(whale)Eistla 1.000  1.000  2.812 0.093580 .  
## s(time11):as.factor(whale)Balder 1.000  1.000 17.759 2.56e-05 ***
## s(time11):as.factor(whale)Jonas  1.001  1.002  0.159 0.690277    
## s(time11):as.factor(whale)Mutti  1.000  1.000  0.369 0.543610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0464   Deviance explained = 6.19%
## UBRE = 1.8611  Scale est. = 1         n = 1203

2.2.2 AIC

AIC(gamwhales1)
## [1] 6732.858

2.2.3 Predicting influence of the different variables

dd1<-sampledata%>%
  summarise(depth=mean(depth,na.rm=T))# creating a data set to predict the influence of time


dd2<-sampledata%>%
  summarise(time11=mean(time11,na.rm=T))# creating a data set to predict the influence of depth

originaldata011<-data.frame(depth=sampledata$depth,time11=mean(sampledata$time11,na.rm=T))# creating a data set to predict the influence of depth


originaldata111<-data.frame(time11=sampledata$time11,depth=mean(sampledata$depth,na.rm=T))# creating a data set to predict the influence of time



datat011<-data.frame(pclick=predict(gamwhales11,newdata = originaldata011,type="response"),depth=originaldata011$depth)# predicting the influence of depth in the number of clicks produced

confintervals<-data.frame(confintervals=predict(gamwhales11,newdata=(originaldata011),type="response",se.fit=T))# get the corresponding CI

datat011$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat011$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit

datat011$variable<-"Depth (m)"# which variable corresponds

datat111<-data.frame(pclick=predict(gamwhales11,newdata = originaldata111,type="response"),time11=originaldata111$time11)# predicting the influence of time in the number of clicks produced

confintervals1<-data.frame(confintervals=predict(gamwhales11,newdata=(originaldata111),type="response",se.fit=T))# get the corresponding CI

datat111$bigfit<-confintervals1$confintervals.fit+confintervals1$confintervals.se.fit
datat111$smallfit<-confintervals1$confintervals.fit-confintervals1$confintervals.se.fit

datat111$variable<-"Time"# which variable corresponds


datat01<-full_join(datat011,datat111)#Joinning both infos about the variables 
## Joining, by = c("pclick", "bigfit", "smallfit", "variable")
datat01$value<-ifelse(datat01$variable=="Depth (m)",datat01$depth,datat01$time11)# turn the variable depth into the real depth

datat01$value<-ifelse(datat01$variable=="Time",datat01$value/(60*60*24),datat01$value)# turn the variable time into the real time in days

datat01$variable1<-ifelse(datat01$variable=="Time","Time since regular echolocation clicks were resumed (days)",datat01$variable)# new labels

2.2.4 Plotting model results/ Figure 5

ggplot(data=datat01,aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+facet_wrap(~variable1,scales="free_x", labeller = labeller(variable1 = label_wrap_gen(32)))+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

p11<-ggplot(data=subset(datat01,variable1=="Time since regular echolocation clicks were resumed (days)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("")+xlab("Time since regular echolocation clicks \nwere resumed (days)")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))+ylim(2,8)

p21<-ggplot(data=subset(datat01,variable=="Depth (m)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))+ylim(2,8)

ggpubr::ggarrange(p21,p11,ncol=2,common.legend = T,widths = c(0.5,0.45))

2.2.5 Predicting number of clicks for the clicking periods from the original dataset

onedata$prediction11<-predict(gamwhales11,newdata=(onedata),type="response")# get the prediction from the model

2.2.6 Plot with the predicted number of clicks for the clicking periods from the original dataset

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction11),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=18),
        plot.title = element_text(size=21),
        strip.text.x = element_text(size = 21))

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction11),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=18),
        plot.title = element_text(size=21),
        strip.text.x = element_text(size = 21))

2.3 GAM number of clicks per second as a function of depth by whale

gamwhales12<-gam(clicks~s(depth,k=7,by=as.factor(whale)),data=sampledata,family=poisson(link="log"))

2.3.1 Summary of the model

summary(gamwhales12)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## clicks ~ s(depth, k = 7, by = as.factor(whale))
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.41232    0.01765   80.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf Ref.df Chi.sq  p-value    
## s(depth):as.factor(whale)Freya  5.421  5.814 18.145 0.003678 ** 
## s(depth):as.factor(whale)Thora  5.764  5.973 25.529 0.000182 ***
## s(depth):as.factor(whale)Mára   2.317  2.867  4.066 0.331349    
## s(depth):as.factor(whale)Frida  2.428  2.959 19.767 0.000158 ***
## s(depth):as.factor(whale)Eistla 2.149  2.683 12.328 0.004252 ** 
## s(depth):as.factor(whale)Balder 4.327  5.078 31.131 1.07e-05 ***
## s(depth):as.factor(whale)Jonas  4.850  5.523 28.508 9.85e-05 ***
## s(depth):as.factor(whale)Mutti  1.000  1.000 19.633 9.96e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0335   Deviance explained = 4.71%
## UBRE = 1.8933  Scale est. = 1         n = 1203

2.3.2 AIC

AIC(gamwhales12)
## [1] 6771.621

2.3.3 Predicting influence of the different expl. variables

originaldata112<-data.frame(depth=sampledata$depth,whale=sampledata$whale)# creating a data set to predict the influence of depth

datat12<-data.frame(pclick=predict(gamwhales12,newdata = originaldata112,type="response"),depth=originaldata112$depth,whale=originaldata112$whale)# predicting the influence of depth in the number of clicks produced

confintervals<-data.frame(confintervals=predict(gamwhales12,newdata=(originaldata112),type="response",se.fit=T))# get the corresponding CI

datat12$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat12$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit

datat12$variable<-"Depth (m)"# which variable corresponds


datat2<-(datat12)#Joinning both infos about the variables 

datat2$value<-ifelse(datat2$variable=="Depth (m)",datat2$depth,datat2$time11)# turn the variable depth into the real depth

datat2$whale<-factor(datat2$whale,levels=c("Freya","Thora","Mára","Frida","Eistla","Balder","Jonas","Mutti"))# correct order of the whale  id

2.3.4 Plotting model results

ggplot(data=subset(datat2,variable=="Depth (m)"),aes(x=value,y=pclick,color=whale,fill=whale))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

2.3.5 Predicting number of clicks for the clicking periods from the original dataset

onedata$prediction12<-predict(gamwhales12,newdata=(onedata),type="response")# get the prediction from the model

2.3.6 Plot with the predicted number of clicks for the clicking periods from the original dataset

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction12),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 17), 
        axis.text.x = element_text(size = 18,angle=60, hjust=1),
        axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 17),
        legend.title = element_text(size=17),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction12),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 8,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 18,angle=60, hjust=1),
        axis.text.y = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=17),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

2.4 GAM number of clicks per second as a function of depth

gamwhales13<-gam(clicks~s(depth,k=7),data=sampledata,family=poisson(link="log"))

2.4.1 Summary of the model

summary(gamwhales13)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## clicks ~ s(depth, k = 7)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.44647    0.01405     103   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##            edf Ref.df Chi.sq p-value    
## s(depth) 3.667  4.466  78.08  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0237   Deviance explained = 2.33%
## UBRE = 1.9234  Scale est. = 1         n = 1203

2.4.2 AIC

AIC(gamwhales13)
## [1] 6807.83

2.4.3 Predicting influence of the different expl. variables

originaldata113<-data.frame(depth=sampledata$depth)# creating a data set to predict the influence of depth

datat13<-data.frame(pclick=predict(gamwhales13,newdata = originaldata113,type="response"),depth=originaldata113$depth)# predicting the influence of depth in the number of clicks produced

confintervals<-data.frame(confintervals=predict(gamwhales13,newdata=(originaldata113),type="response",se.fit=T))# get the corresponding CI

datat13$bigfit<-confintervals$confintervals.fit+confintervals$confintervals.se.fit
datat13$smallfit<-confintervals$confintervals.fit-confintervals$confintervals.se.fit

datat13$variable<-"Depth (m)"# which variable corresponds


datat3<-(datat13)#Joinning both infos about the variables 

datat3$value<-ifelse(datat3$variable=="Depth (m)",datat3$depth,datat3$time11)# turn the variable depth into the real depth

2.4.4 Plotting model results

ggplot(data=subset(datat3,variable=="Depth (m)"),aes(x=value,y=pclick))+geom_line(size=1.2)+ylab("Predicted number of clicks / s \nfor the clicking periods")+xlab("Depth (m) \n")+
  geom_ribbon(aes(ymin = smallfit,ymax=bigfit,x=value), alpha = .4)+
  theme_bw()+  scale_color_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+  scale_fill_manual(values=c("black", "blue4","dodgerblue1","turquoise4","green4","darkorchid4","coral4","darkorange"))+labs(fill="",color="")+
  theme(
        axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 19),
        axis.text.y = element_text(size = 19),
        axis.title.y = element_text(size = 20),
        legend.title = element_text(size=20),
        legend.text= element_text(size=19),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

2.4.5 Predicting number of clicks for the clicking periods from the original dataset

onedata$prediction13<-predict(gamwhales13,newdata=(onedata),type="response")# get the prediction from the model

2.4.6 Plot with the predicted number of clicks for the clicking periods from the original dataset

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction13),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

ggplot(originaldata, aes(y=-Depth, x=sst/(60*60*24)))+
    geom_point(color="black",shape=".") +
   geom_point(data=onedata,aes(colour=(prediction13),y=-depth,x=sst/(60*60*24)),alpha=0.3,shape=".")+
  xlab("Time since regular echolocation clicks were resumed (days)")+ylab("Depth(m)")+facet_wrap(~whale,scales="free_x",ncol=4,nrow=2)+labs(color="Predicted number \nof clicks \\ s")+scale_color_gradient2(low = "blue",
midpoint = 5.5,
mid = "green",
high = "red",
space="Lab",
  na.value = "black")+
  theme(
        axis.title.x = element_text(size = 18), 
        axis.text.x = element_text(size = 17,angle=60, hjust=1),
        axis.text.y = element_text(size = 17),
        axis.title.y = element_text(size = 18),
        legend.title = element_text(size=18),
        legend.text= element_text(size=17),
        plot.title = element_text(size=20),
        strip.text.x = element_text(size = 21))

2.5 Comparing models

data.frame("Response variable"=rep("Number of clicks per second",4),"Explanatory variables"=c("depth by narwhal ID and time by narwhal","depth and time","depth by narwhal ID","depth"),"AIC"=c(round(AIC(gamwhales1),3),round(AIC(gamwhales11),3),round(AIC(gamwhales12),3),round(AIC(gamwhales13),3)))

2.5.1 Estimating the cue rate using the model results

tagdur<-originaldata%>%
  dplyr::group_by(whale)%>%
  dplyr::summarise(sumone=unique(len_click),ll=unique(tot_len))


onedata<-full_join(onedata,tagdur)
## Joining, by = "whale"
est_cr<-onedata%>%
  group_by(whale)%>%
  summarise(meany=mean(prediction1,na.rm=T),sumone=unique(sumone),ll=unique(ll),mean=(meany*sumone)/ll,sd=sd(prediction1,na.rm=T),cv=sd/meany,
              lower = meany-qt(0.975,df=n()-1)*sd/sqrt(n()),
              upper = meany+qt(0.975,df=n()-1)*sd/sqrt(n()))

2.5.1.1 Results

n <- as.numeric(nrow(est_cr))
xbar <- mean(est_cr$mean,na.rm=T) 
s <- sd(est_cr$mean,na.rm=T)

#calculate margin of error
margin <- qt(0.975,df=n-1)*s/sqrt(n)

#calculate lower and upper bounds of confidence interval
low <- xbar - margin
high <- xbar + margin


cv0<-s/xbar
  
data0<-data.frame(mean=xbar,low=low,high=high,cv=cv0,sd=s,id="Standard Average")
#var(df5$mean)

The standard average for the cue rate is 1.412 and the corresponding 95% CI is [0.998 ,1.827]. The cv obtained is 0.351 and the sd is 0.495.

2.5.1.2 Weighted Results

wxbar<-  weighted.mean(est_cr$mean,est_cr$ll,na.rm=T) 
ws <- sqrt(var.wtd.mean.cochran(est_cr$mean,est_cr$ll))

#calculate margin of error
wmargin <- qt(0.975,df=n-1)*ws

#calculate lower and upper bounds of confidence interval
wlow <- wxbar - wmargin
whigh <- wxbar + wmargin

wcv<-ws/wxbar
#var(df5$mean)

wdata0<-data.frame(mean=wxbar,low=wlow,high=whigh,cv=wcv,sd=ws,id="Weighted Average")

The weighted average for the cue rate is 1.29 and the corresponding 95% CI is [0.961 ,1.62]. The cv obtained is 0.108 and the weighted sd is 0.139.

2.5.1.3 Joining results from cue rates

adata<-full_join(data0,wdata0)
## Joining, by = c("mean", "low", "high", "cv", "sd", "id")
adata$pred<-"Prediction"
data$pred<-"Original"
cr_data<-full_join(adata,data)
## Joining, by = c("mean", "low", "high", "cv", "sd", "id", "pred")
cr_data<-data.frame(cr_data[cr_data$id%in%c("Weighted Average","Standard Average"),])
data.frame(cr_data[cr_data$id%in%c("Weighted Average","Standard Average"),c("id","pred","mean","sd","cv","low","high")])

2.5.1.4 Plotting results

ggplot(cr_data, aes(x= id, y= mean,color=pred,fill=pred))+
 # geom_point()+ 
  geom_pointrange(aes(ymin=low,ymax=high,color=pred,fill=pred), 
  position = position_dodge(width = 0.5))+
  ylab("Mean Cue Rate")+
  xlab("")+
  labs(color="",fill="")+
  theme(axis.text.x=element_text(angle=30, hjust=1))

2.5.2 Saving the best model and its results

save(file="model_results_n_clicks.rda",list=c("gamwhales1","datat"))